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ABSTRACT 

We present a new technique for calibrating the primary beam of a wide-field, drift-scanning antenna 
element. Drift-scan observing is not compatible with standard beam calibration routines, and the 
situation is further complicated by difficult-to-parametrize beam shapes and, at low frequencies, the 
sparsity of accurate source spectra to use as calibrators. We overcome these challenges by building up 
an interrelated network of source "crossing points" — locations where the primary beam is sampled 
by multiple sources. Using the single assumption that a beam has 180° rotational symmetry, we can 
achieve significant beam coverage with only a few tens of sources. The resulting network of crossing 
points allows us to solve for both a beam model and source flux densities referenced to a single 
calibrator source, circumventing the need for a large sample of well-characterized calibrators. We 
illustrate the method with actual and simulated observations from the Precision Array for Probing 
the Epoch of Reionization (PAPER) . 

Subject headings: instrumentation: interferometers — methods: miscellaneous — techniques: inter- 
ferometric 



1. INTRODUCTION 

The past decade has seen a renewed interest in low fre- 
quency radio astronomy with a strong focus on cosmol- 
ogy with the highly redshifted 21cm line of neutral hy- 
drogen. Numerous facilities and experiments are already 
online or under construction, includi ng the Gia nt Metre- 
Wave Radio Telescope (GMRT; Swarup et all 1991 f\ the 
LOw Frequency ARray (LOFAR; Rottgcring 20030 the 
Long Wavelength Array (LWA; |Ellingson et al.||2009| p1 
and the associated Large Aperture-experiment to Detect 
the Dark Ages (LEDA) experiment, the Cylindri cal Ra- 



son et al. 2006 Seo et al. 2010 



dio Telescope (CRT/BAORadio, formerly HSHS, Peter- 



the Experiment to 



Detec t the Global EoR Step (EDGES; Bowman et al 



2008 



t he Mu rchison Widefield Array (MWA; Lons 



dale et al. 
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and the Donald C. Backer Precision 
Array toT^Probing the Epoch of Reionization (PAPER; 
Parsons et al. [20101}^ 21cm cosmology experiments will 
need to separate bright galactic and extragalactic fore- 
grounds from the neutral hydrogen signal, which can be 
faint er by as much as 5 orders of magnitude or more (see, 
e.g., Furlanetto et al. 2006, and |Santos et al. 2005). As 
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such, an unprecedented level of instrumental calibration 
will be necessary for the detection and characterization 
of the 21cm signal. 

Achieving this level of calibration accuracy is compli- 
cated by the design choice of many experiments to em- 
ploy non-tracking antenna elements (e.g. LWA, MWA, 
LOFAR and PAPER). Non-tracking elements can pro- 
vide significant reductions in cost compared to tradi- 
tional dishes, while also offering increased system stabil- 
ity and smooth beam responses. However, non-tracking 
elements also present many calibration challenges be- 
yond those of traditional radio telescope dishes. Most 
prominently, the usual approach towards primary beam 
calibration — pointing at and dithering across a well- 
characterized calibrator source — is not possible. In- 
stead, each calibrator can only be used to characterize 
the small portion of the primary beam it traces out as 
it passes overhead. Additionally, the wide fields of view 
of many elements make it non-trivial to extract individ- 
ual cal ibrator sources f rom the data ( see, e.g., Parsons 
fc Backer 2009 and Paciga et al.||2011| for approaches to 



isolate calibrator sources) . Finally, many of these arrays 
use dipole and tile elements, the response of which are 
not easily described by simple analytic functions. 

In this paper, we present a method for calibrating the 
primary beams of non-tracking, wide-field antenna ele- 
ments using astronomical sources. We illustrate the tech- 
nique using both simulated and observed PAPER data 
from a 12 antenna array at the NRAO site in Green 
Bank, WV. PAPER is an interferometer operating be- 
tween 100 and 200 MHz, targeted towards the highly- 
redshifted 21cm signal from the epoch of reionization. 
Although the cosmological signal comes from every di- 
rection on the sky, an accurate primary beam model will 
be necessary to separate the faint signal from bright fore- 
ground emission. In this work, we use a subset of the 
brightest extragalactic sources to calibrate the primary 
beam. Because interferometers like PAPER are insensi- 



2 



tive to smooth emission on large scales, these extragalac- 
tic sources are easily detectable despite the strongly in- 
creasing brightness of Galactic synchrotron emission at 
low radio frequencies. Only the measured relative flux 
densities of each source are needed to create a beam 
model, but to facilitate comparision with other catalogs, 
we us e the a bsolute spectrum of Cygnus A from |Baars| 
et al. (1977) to place our source measurements on an 
absolute scale. 

The structure of this paper is as follows: in f|2]we mo- 
tivate the problem and the need for a new approach to 
primary beam calibration for wide-field, drift-scanning 
elements. In ^|3] we present our technique for primary 
beam calibration. We show the results of applying the 
method to simulated and actual observations in and 
^ respectively, and we conclude in ^ 

2. MOTIVATION 

For non-tracking arrays with static pointings such as 
PAPER, every celestial source traces out a repeated 
"source track" across the sky, and across the beam, each 
sidereal day. The basic relationship between the per- 
ceived source flux density (which we shall call a measure- 
ment, M) measured at time t and the intrinsic source flux 
density (/) is: 

M{t) = b{s{t))f, (1) 

where b is the response of the primary beam toward the 
time-dependent source location s. 

If the inherent flux density of each source were well- 
known, it would be straightforward to divide each M{t) 
by / to obtain b along the source track s{t). To form 
a complete beam model, one would then need enough 
well-characterized sources to cover the entire sky. In 
the 100-200 MHz band, however, catalog accuracy for 
most sources is lagging behind the need f or precise beam 
calibration (Vollnicr ct al.^.2G05j ^Jacobs et al. 2011), 
which in turn is necessary tor generating improved cata- 
logs. Without accurate source flux densities, both b and 
/ are unknowns, and Equation [T] is underconstrained. 
The problem becomes tractable if several sources pass 
through the same location on the sky, and therefore are 
attenuated by the same primary beam response. How- 
ever, the density of bright sources at low frequencies is 
insufficient to relate sources at different declinations. Ad- 
ditional information is necessary to break the degeneracy 
between primary beam response and the inherent flux 
densities of sources. 

3. METHODS 

One way to overcome the beam-response/flux-density 
degeneracy described above is to assume 180° rotational 
symmetry in the beam. Under this assumption, each 
source creates two source tracks across the sky: one cor- 
responding to the actual position of the source, and the 
other mirrored across the beam center. Under this sym- 
metry assumption, tracks overlap at "crossing points," 
as schematically illustrated in Figure [l] for two sources: 
Cygnus A and Virgo A. At a crossing point, there are 
only 3 unknowns, since each source is illuminated by the 
same primary beam response, b. Since two circles on the 
sky cross twice (if they cross at all), there will be two 
independent relations that together provide enough in- 
formation to constrain both the primary beam response 




Figure 1. The building up of source traclss and crossing points. 
Left: Tlie paths Cygnus A (top) and Virgo A (bottom) take across 
the beam for the PAPER array in Green Bank. The projection 
is orthographic with zenith at the center; dotted lines are 10° and 
30° steps in elevation and azimuth, respectively. Right: The tracks 
of Cygnus A and Virgo A, and their 180° rotated images overlaid. 
There are now 6 crossing points. These are the locations where the 
beam and source parameters can be solved for independently. 

at those points and the flux densities of sources passing 
through them (relative to an absolute flux scale). By 
using multiple sources, it is possible to build a network 
of crossing points that covers a large fraction of the pri- 
mary beam. Furthermore, such a network allows data 
to be calibrated to one well-measured fiducial calibra- 
tor. For observations with the Green Bank deployment 
of PAPER, source fluxes are related to Cygnus A from 



Baars et al.| ( |1977[ ). 

In the rest ot this section, we discuss the algorithmic 
details of how this approach to primary beam calibration 
is implemented. These steps are: 

1. extracting measurements of perceived so urce flux 
densities versus time from observations ([3.1), 



2. entering measureme nts in to a gridded sky and find- 
ing crossing points (^3.2[), 



3. solving a least -squares matrix formulation of the 
problem (^3.3), and, finally, 



4. deconvolving the irreg ularly sampled sky to create 
a smooth beam (^3.4). 



We also discuss prior inform ation that can be included 
to refine the beam model in ^3.5 



3.1. Obtaining Perceived Source Flux Densities 

The principal data needed for this approach are mea- 
surements of perceived flux density versus time for mul- 
tiple calibrator sources as they drift through the primary 
beam. In practice, any method of extracting perceived 
individual source flux densities (such as image plane anal- 
ysis) could be used; the beam calibration procedure is 
agnostic as to the source of these measurements. In this 
work, we use delay/delay-rate (DDR) filters to extract 
estimates of individual source fiux densities as a func- 
tion of time and frequency. The frequency information 
can be used to perform a frequency-dependent beam cal- 
ibration if the SNR in the observations is high enough. 
These filters work in delay /delay-rate space (the Fourier 
dual of frequency/time space) to isolate fiux density per 
baseline from a specific geometric area on the sky; for a 
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full description of the tec hnique, the reader is referred to 
Parsons fc Backer] ( [2009] ) . 

'i'he hrst step of our pipeline to produce perceived flux 
density estimates is to filter the Sun from each baseline 
individually in DDR space. This is done to ensure that 
little to no signal from the Sun, which is a partially- 
resolved and time-variable source, remains in the data. 
(In principle, data from different observing runs sepa- 
rated by several months could provide complete sky cov- 
erage while avoiding daytime data altogether. We chose 
to use data from only one 24 hour period to minimize the 
effect of any long timescale instabilities in the system). 
A Markov-Chain Monte Carlo (MCMC) method for ex- 
tracting accurate time- and frequency-dependent source 
models via self-calibration using DDR-based estimators 
is then used to model the 4 brightest sources remaining: 
Cygnus A, Cassiopeia A, Virgo A, and the Crab Nebula 
(Taurus A). The MCMC aspect of this algorithm iter- 
ates on a simultaneous model of all sources being fit for 
to minimize sidelobe cross-terms between sources. After 
removing the models of Cygnus A, Cassiopeia A, Virgo 
A and Taurus A from the data, a second pass of the 
MCMC DDR self-calibration algorithm extracts models 
of the remaining 22 sources listed in Table [T] 

3.2. Gridding the Measurements 

We can increase the signal-to-noise at a crossing point 
by combining all measurements within a region over 
which the primary beam response can be assumed to be 
constant. To define these regions, we grid the sky. In this 
work, the beam model i s constructed on an HEALPix 
map ( [Gorski et al. 2005) with pixels ~ 0.9° on a side. 
The choice ot grid pixel size is somewhat arbitrary. Us- 
ing a larger pixel size broadens the beam coverage of each 
source track, creating more crossing points and helping 
to constrain the overall beam shape. However, when the 
pixels are too large, each pixel includes data from sources 
with larger separations on the sky. Since the principal 
tenet of this approach is that each source within a cross- 
ing point sees the same primary beam response, exces- 
sively large pixels can violate this assumption, result- 
ing in an inaccurate beam model. For PAPER data, a 
HEALPix grid with 0.9° pixels is found to be a good 
balance bet ween these competing factors, as will be ex- 
plained in j |4.1[ For other experiments with narrower, 



more rapidly evolving primary beams, smaller pixels may 
be necessary. 

To introduce our measurements of perceived flux den- 
sity into the grid, we first recast Equation [l] into a dis- 
crete form: 

= bjk, (2) 

where Mi and bi are the respective perceived source flux 
densities and primary beam responses in the pixel i, and 
fc is a source index labelling the inherent source flux den- 
sity, /. To generate a single measurement of a source 
for each pixel i, we use a weighted average of all mea- 
surements of that source falling with a single pixel. The 
weights are purely geometric and come from interpolat- 
ing the measurement between the four nearest pixels. 

3.3. Forming a Least-Squares Problem 

Once all the data are gridded, we solve Equation [2] for 
all crossing pixels simultaneously. To do this, we set up 
a linearized least-squares problem using logarithms: 



log(M,) = log(6,) + log(/fe). 



(3) 



Because thermal noise in measurements becomes non- 
Gaussian in Equation [3j the solution to the logarithmic 
form ulation of the least-squares problem is biased. In 
§4.I[ we investigate the effect of this bias using simu- 
lations and find that the accuracy of our results is not 
limited by this bias, but instead by sidelobes of other 
sources. Therefore, while Equation [T] can in principle be 
solved without resorting to logarithms using an iterative 
least-squares approach, we find that this is not necessary 
given our other systematics. 

Once the logarithms are taken, we can construct a solv- 
able matrix equation, which can be expressed generally 
as: 

WM = WAX, (4) 

where W is a column- matrix of weights, M is a column- 
matrix of measurements, A is the matrix-of-condition, 
describing which measurements are being used to con- 
strain which parameters, and X is a column matrix con- 
taining all the parameters we wish to measure: the beam 
responses b and the source flux densities /. To illustrate 
the form of this equation, we present the matrix rep- 
resentation of the Cygnus A/Virgo A system shown in 
Figure [TJ 
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On the left-hand side of the equation are the two col- 



umn matrices W and M. The weights, W^, are defined 
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in §3.3. 1| M contains logarithms of the perceived source 
flux density measurements in each pixel, Mi. Recall that 
each Mi corresponds to one source. 

On the right-hand side, the weighting column matrix 
W appears again, followed by the matrix-of-condition 
A, and then X, a column matrix containing the pa- 
rameters we wish to solve for: the primary beam re- 
sponse at the 6 crossing points and the flux densities of 
Virgo and Cygnus. The matrix-of-condition, A, iden- 
tifies which sources and crossing points are relevant for 
each equation. We have schematically divided it: to the 
left of the vertical line are the indices used for select- 
ing a particular crossing point; to the right are those for 
the sources. The first 4 lines represent the two northern 
Cygnus/Virgo crossing points, and the next four repre- 
sent the two southern ones (which are identical copies 
of the northern ones) . Finally, the last 4 lines represent 
the 2 points where Virgo crosses itself; notice that the 
Cygnus source column is blank for these 4 rows. 

It should be noted that Equations [4] and [5] contains no 
absolute flux density reference. The simplest way to set 
this scale is to treat all the recovered flux densities as 
relative values compared to the flux calibrator (Cygnus 
A in the case presented here) . One can then place all the 
flux densities onto this absolute scale. An equally valid 
approach is to append an extra equation with a very high 
weight, which sets the flux calibrator to its catalog flux 
value. 

3.3.1. Weighting of Measurements m the Least-Squares 
Formalism 

In a least-squares approach, optimal weights are in- 
versely proportional to the variance of each measure- 
ment. To calculate the variance of each measurement, 
we must propagate the uncertainty in the initial per- 
ceived flux density measurements through the averaging 
and logarithm steps. The noise level in each interfero- 
metric visibility is roughly constant, and the DDR filters 
average over a fixed number of visibilities for each per- 
ceived flux density estimate, leading to equal variance at 
each time sample. To produce optimal weights account- 
ing for the many time samples averaged into each beam 
pixel and the propagation of noise through the logarithm, 
each logarithmic measurement in Equation [4] should be 
weighted by: 



(6) 



where j indexes the time step (which is generally fast 
enough to produce many measurements inside a pixel), 
Wj is the geometric sub-pixel interpolation weighting of 
each measurement, and Mi is the weighted average of 
all measurements of a source's fiux density in pixel i. 
Without the square root, these weights would be propor- 
tional to the inverse of the variance in each logarithmic 
measurement; the sum over the geometric weights is the 
standard reduction of variance for a weighted average, 
and the factor of Mf comes from propagating variances 
through the logarithm. The square root appears because 
a least-squares solver using matrix inversion will add in 
an additional factor of the weight, leading to the desired 
inverse-variance weighting. Other solvers using different 
methods may require different weights. 




Figure 2. The sky coverage of crossing points and source tracks. 
Left: The location of crossing points for the 25 sources listed in 
Table [T] used to calibrate the PAPER primary beam. The projec- 
tion is orthographic with zenith at the center; dotted lines are 10° 
and 30° steps in elevation and azimuth, respectively. Right: The 
sky coverage of all 25 source tracks. Although the least-squares 
inversion solves only for the beam response at crossing points, we 
can include all source track data using the recovered flux density 
of each source to create a primary beam estimate along the entire 
track. 

3.3.2. Solving the Least-Squares Equation for a First-Order 

Beam Model 

The solution of Equation|4j the matrix A, contains two 
distinct sets of parameters: the beam responses at cross- 
ing points and the flux densities of each source (mod- 
ulo an absolute scale). There is additional information 
that may be included to improve the beam model beyond 
that generated by solving for the responses at crossing 
points. Given the flux-density solutions for each source, 
each source track now provides constraints on beam pix- 
els that are not crossing points. By dividing each per- 
ceived flux density source track by the estimated inherent 
flux density from the least-squares inversion, we produce 
a track of primary beam responses, with greater cover- 
age than that provided by crossing points alone. We 
illustrate the difference in coverage in Figure [2j The left 
hand panel shows the locations of crossing points for the 
25 sources used to calibrate the PAPER beam, listed in 
Table [T] The right hand panel shows the increased beam 
coverage that comes from including non-crossing point 
source-track data. To create an initial beam model, we 
average within each pixel the estimated beam responses 
from each source, weighting by the estimated flux of that 
source. Given equal variance in each initial perceived flux 
density measurement, this choice of weights is weighting 
by signal to noise. 

3.4. Using Deconvolution to Fill in Gaps in the Beam 

Model 

To produce a model of the beam response in any direc- 
tion, we must fill in the gaps left by lim ited sky sampling. 
We use a CLEAN-like ( |H6gbom|[l974l ) deconvolution al- 
gorithm in Fourier space to fill in the holes in the beam. 
We iteratively fit and remove a fraction of the bright- 
est Fourier component of our beam until the residuals 
between the model and the data are below a specified 
threshold. In addition to measured beam responses de- 
rived from source tracks, we add the constraint that the 
beam must go to zero beyond the horizon. In the decon- 
volution, each pixel is weighted by the estimate of the 
beam response in that pixel, reflecting that SNR will be 
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highest where the beam response is largest. This beam- 
response weighting was unnecessary in previous steps, 
since the least-squares approach solved for each pixel 
independently. This weighting scheme again represents 
signal-to-noise weighting, given the equal variance in our 
initial perceived flux density estimates. The result of 
the deconvolution is an interpolated primary beam model 
with complete coverage across the sky. 

3.5. Introduction of Prior Knowledge 

Up to this point, we have made only two fairly weak 
assumptions about our beam: that it possesses 180° rota- 
tional symmetry, and that the response is zero below the 
horizon. However, if we have additional prior knowledge 
about our beam, we can better constrain the final model. 
In particular, we can use beam smoothness constraints 
to identify unphysically small scale features introduced 
by sidelobes in the source extraction or by incomplete 
sampling in the deconvolution. We choose to incorpo- 
rate additional smoothness information by filtering our 
model in Fourier space to favor large-scale modes. 

PAPER dipoles were designed with emphasis on spa- 
tial and spectral smoothness in primary beam shape. 
To smooth the PAPER beam model, we choose a cutoff 
in Fourier space that corresponds to the scale at which 
> 99.9% of the power is accounted for in a computed 
electromagnetic model of the beam. While such a filter is 
not necessarily generalizable to other antenna elements, 
we find it necessary to suppress the substantial sidelobes 
associated with observations from a 12-antenna PAPER 
array that are discussed below. 

4. APPLICATION TO SIMULATED DATA 

To test the robustness of this approach, we apply it to 
several simulated data sets. The results of these simu- 
lations, with and without Gaussian noise, are described 
below in |4.1| We also simulate raw visibility data to test 
the effectiveness of s ource extraction; these simulations 
are discussed in S 4.2 The major difference between these 



two methods of simulation is that the simulation using 
raw visibility data allows for imperfect source isolation, 
leading to contamination of the source tracks by side- 
lobes of other sources. These sidelobes have a significant 
effect on the final beam model that is derived. 

4.1. Simulations of Perceived Flux Density Tracks 

We simulate perceived flux density tracks using several 
model beams of different shapes, including ones with sub- 
stantial ellipticity and an ~ 15° rotation around zenith. 
We also input several source catalogs, including a case 
with 10, 000 Jy sources spaced every degree in declina- 
tion, and a case using the catalog values and sources 
listed in Table [T] that approximately match the sources 
extracted from observations. In all combinations of beam 
models and source catalogs, we recover the input beam 
and source values with < 5% error. The average error in 
source flux density is 2.5%. We see no evidence for resid- 
ual bias, as the distribution of error is consistent with 
zero mean to within one standard deviation. 

We also test the effect of adding various levels of Gaus- 
sian noise to a simulation involving a fiducial beam model 
and the 25 sources listed in Table [l] Only when the noise 
level exceeds an RMS of 10 Jy in each perceived flux den- 
sity measurement does the mean error in the solutions 



exceed 10%. As the noise increases beyond this level, 
there is a general trend to bias recovered flux densities 
upward; as mentioned earlier, this bias is introduced by 
the logarithms in EquationjSj However, the expected cor- 
responding noise level of DDR-extracted perceived flux 
density measurements from our 12-element PAPER ar- 
ray is < 1 Jy per sample. At a simulated rms noise level 
of 1 Jy, the solutions are recovered with a mean error of 
< 3%. This result validates our previous statement that 
the bias introduced by the logarithms in Equation |3] is 
not a dominant source of error. 

It is also worth noting that these simulations were used 
to identify 0.9° as the best HEALPix pixel side for our 
grid; with too large a pixel size (1.8°) the model becomes 
significantly compromised. This results from combining 
measurements of sources that are subject to significantly 
different primary beam responses. We choose not to use 
a smaller pixel size, since it will reduce the fractional sky 
coverage of our crossing points and increase the compu- 
tational demand of the algorithm. 

4.2. Simulations of Visibilities 

We also apply this technique to simulated visibilities 
in order to test the complete analysis pipeline, includ- 
ing the DDR-based estimation of perceived source flux 
densities. The visibility simulations are implemented in 
the AIP^vF^ software toolkit. The simulated observa- 
tions correspond to actual observations made with the 
12-element PAPER deployment in Green Bank described 
in Sj5] We match the observations in time, antenna posi- 
tion, and bandwidth. We also include the expected level 
of thermal noise in the simulated visibilities. We simulate 
"perceived" visibilities by attenuating the flux density of 
each source by a model primary beam. The DDR filters 
return estimates of perceived source flux density for the 
input primary beam. 

In these simulations, we only include bright point 
sources and a uniform-disk model of the Sun. As a re- 
sult, source extraction is expected to be more accurate in 
simulation than in real data, since the sources we extract 
account for 100% of the simulated signal. However, these 
simulations do provide a useful test of DDR-fllter-based 
source extraction and of the level of contamination from 
sidelobes. 

As might be expected, the estimates of perceived 
source flux density versus time from simulated visibilities 
contain structure that is not attributable to the beam. 
These features are almost all due to sidelobes of other 
sources; they persist in real data and are reproducible 
over many days. Figure [3] shows the difference in the 
source tracks produced for Cassiopeia A by each simu- 
lation method (cases 1 and 2 in the flgure), as well as 
the track extracted from the observed data described ^}5] 
(case 3). We find that the DDR filter only extracts a 
fraction of the flux density from the Crab Nebula. This 
bias is unique to Crab, and is most likely a result of the 
proximity of the Sun to Crab in these observations, cou- 
pled with the limited number of independent baselines 
in the 12-antenna array. For this reason, we choose to 
exclude the Crab from our final analysis. 

Pe rforming the least-squares inversion described in 
93.31 on the simulated source tracks, we find that the 
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Cas A Perceived Flux Tracks 



Case 1 (no sidelobes) 
Case 2 {with sidelobes) 
Case 3 (observed data) 




5 10 15 20 25 

Local Sidereal Time (Hours) 

Figure 3. Three different perceived flux tracks for Cassiopeia A. 
Case 1 is simply a cut through a model beam multiplied by a cata- 
log flux density value; this is the simplest case and contains no side- 
lobe contamination. Case 2 uses DDR filters to extract the source 
track from simulated visibilities; case 3 does the same, but applied 
to observed data. The simulated visibilities in case 2 demonstrate 
the features that arise from sidelobes of other sources during the 
DDR source extraction. However, the simulation employs a sim- 
plistic point-source sky and cannot reproduce all sidelobe features 
seen in the real data. 




19 20 21 

Local Sidereal Time (Hours) 

Figure 4. The effect of the Fourier domain filter. The deconvolu- 
tion of i ]3.4| fllls in gaps in sky coverage, but leaves small scale struc- 
ture that IS not associated with the primary beam (gray curve). 
By flltering the solution in the Fourier domain to incorporate prior 
knowledge of beam smoothness, one can achieve an improved beam 
model (black curve). 



sidelobes introduce errors into estimates of source flux 
densities. The average error in flux density is 10%, with 
the largest errors exceeding 20%. The distribution of er- 
rors is consistent with zero mean, indicating no strong 
biasing of the flux densities. 

We also find that the resultant primary beam model 
matches the input beam to within 15% percent. How- 
ever, the model contains small scale variations, as seen 
in Figure [4j Using prior knowledge of beam smoothness, 
we can reject these features as unphysical and apply a 

This filter re- 



Fourier-domain filter, as described in ^3.5 
duces the effect of sidelobes in the source tracks, as they 
appear on Fourier scales that are not allowed by smooth- 
ness constraints. Even with a relatively weak prior on 
smoothness (i.e., retaining more Fourier modes than are 
present in our input model), we substantially reduce the 
small scale variations in our beam. The effect of the fil- 
ter on small scale variations is illustrated in Figure |4j 
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Figure 5. The 12 element PAPER array in Green Bank. Shaded 
circles represent vertical height relative to antenna 0; the radius of 
the circle is 10 times the height in meters. Light gray is up, dark 
gray is down. 

The filter also reduces large scale errors in the model, 
bringing the overall agreement to with 10% of the input 
beam. 

In summary, we find that although interfering sidelobes 
from other sources compromise our source flux density 
measurements, our approach recovers the input primary 
beam model at the 15% level. With the introduction 
of a Fourier space filter motivated by prior knowledge 
of the beam smoothness, the output model improves to 
within 10% of the input. It is also worth noting that the 
effectiveness of the beam calibration technique presented 
here will substantially improve with larger arrays, where 
sidelobe interference will be reduced and source tracks 



will more closely resemble the ideal case discussed in ^ 4.1 



5. OBSERVED DATA 

In this section, we describe the application of this tech- 
nique to 24 continuous hours of data taken with the 
PAPER array deployed at the NRAO site near Green 
Bank from July 2 to July 3, 2009. At this time, the 
array consisted of 12 crossed-dipole elements. Only the 
north/south linear polarization from each dipole was cor- 
related. The data used in this analysis were observed be- 
tween 123 to 170 MHz (avoiding RFI) and were split in 
420 channels. The array was configured in a ring of 300m 
radius. The longest baselines are 300m, while the short- 
est baseline is 10m; the configuration is shown in Figure 

tThis configuration gives an effective image plane res- 
ution of 0.4°. 

5.1. Data Reduction 

In addition to the DDR-filter algorithm described in 
§3.11 several other pre-processing steps are necessary. 
Per-antenna phase and amplitude calibration are per- 
formed by fringe fitting to Cygnus A and other bright 
calibrator sources. We also use a bandpass calibration 
based on the spectrum of Cygnus A. This calibration is 
assumed to be stable over the full 24 hours. Two further 
steps in the reduction pipeline were first described in 
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278 280 282 284 286 288 290 292 312 314 316 318 320 322 324 
Cable Temperature (K) PDA Temperature (K) 

Figure 6. The effect of ambient temperature on the system gain, 
derived from observations of Cygnus A near zenith, where beam 
effects are unimportant. The left-hand side shows the effect of 
the cable temperature; the right shows pseudo-differential amplifier 
(PDA)/balun temperature. The solid black line is the best fit to 
the data. Different shades of gray represent data from different 
antennas. 



Parsons et al. (2010): gain linearization to mitigate data 
quantization effects in the correlator, and RFI excision. 

New to this work is the use of the cable and balun 
temperatures to remove time-dependent gains caused by 
temperature fluctuations. A thirteenth dipole was op- 
erated as the "gain-o-meter" described in [Parsons et al.| 
(2010). In brief, the "gain-o-meter" is an antenna where 
the balun is terminated on a matched load, rather than a 
dipole. The measured noise power from this load tracks 
gain fluctuations in the system. We record the temper- 
ature of several compo nents of our "gain-o-meter" using 
the system described in Parashare & Bradley (2009). We 
find a strong correlation between the absolute gam of our 
system and the temperature of these components; this ef- 
fect is illustrated in Figure [6] We correct for these gain 
changes by applying a linear correction derived from the 
measured temperature values. This step is crucial for 
the success of the beam calibration; without correcting 
for them, these gain drifts are indistinguishable from an 
east /west asymmetry in the primary beam. After these 
steps, we process the data with the DDR-filtering al- 
gorithm to produce estimates of perceived flux density 
versus time for the 25 sources listed in Table [T] 

5.2. Results 
5.2.1. PAPER Primary Beam Model 

The first estimate of the PAPER primary beam derived 
with this approach is shown in Figure [Tj This figure is 
produced by dividing each source tracK by an estimate 
of its inherent flux density produced by the least squares 
inversion. This transforms each track into an estimate of 
the primary beam response, which are the added together 
into a HEALPix map. There are significant fluctuations 
from pixel to pixel that are clearly unphysical, and relate 
to source sidelobes. There are also significant gaps in 
beam coverage, even with 25 source tracks. 

Figure |8] show the results after deconvolving the sam- 
pling pattern from Figure [7] The beam is now complete 
across the entire sky, although there is still substantial 
small-scale structure unrelated to the inherent beam re- 
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Figure 7. The initial solution for the PAPER primary beam 
response along source tracks. Each perceived fiux density track 
is divided by the estimate of that source's inherent flux density 
produced by the least-squares inversion. This yields estimates of 
the primary beam response along each track. With this initial 
solution, a large fraction of the sky is already covered by the data, 
although there is a fair bit of variation from point to point. The 
color scale is linear and normalized to 1.0 at zenith; as in Figure 
^ dotted lines are 10° and 30° steps in elevation and azimuth, 
respectively. 
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Figure 8. The deconvolved PAPER beam model. The grayscale 
colors are linear and show the primary beam response, normal- 
ized to 1.0 at zenith. Solid black lines are contours of constant 
beam response in 10% increments; the thick black line is the half- 
power point. As in Figure [T] dotted lines are 10° and 30° steps 
in elevation and azimuth, respectively. This model was produced 
by deconvolving the sampling pattern from Figure [7] interpolating 
over the gaps in declination where there are no strong calibrator 
sources. 
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Figure 9. The smooth PAPER beam model. As in Figure [8] the 
grayscale colors are linear and show the primary beam response, 
normalized to 1.0 at zenith. Solid black lines are contours of con- 
stant beam response in 10% increments; the thick black line is the 
half-power point. Dotted lines are 10° and 30° steps in elevation 
and azimuth, respectively. This smooth model was produced by 
retaining only the large scale Fourier components from Figure Is] 
The cutoff in Fourier space was derived from a computed electro- 
magnetic PAPER beam simulation. This Fourier mode selection 
substantially smooths out variations in the initial solution that are 
unphysical based on the smoothness of beam response in the sim- 
ulation. 

sponse. We suppress these features in our final model 
by retaining only the large scale Fourier components, as 
described in §3.4| This final model is shown in Figure 
[9l It is clear that with our deconvolution and choice of 
Fourier modes, we recover a very smooth, slowly evolving 
beam pattern, as expected from previous models of the 
PAPER primary beam. 

5.2.2. Source Catalog 

The measured flux densities of the calibrator sources, 
produced by the least-squares inversion, are presented in 
Table [T] The overall flux scale is normalized to the value 



of Cygnus A reported in Baars et al. (1977]). Error-bars 
are estimated by measuring the change in nux necessary 
to increase Ax^ between the model and the measured 
per ceiv ed source flux density tracks by 1.0. As noted 
in §4.2[ these flux densities can be compromised by side- 
lobes of other sources. However, the results are generally 
accurate to within ±10%. 

5.3. Tests of Validity 

We perform several tests to investigate the validity of 
our results. The comparison of our recovered fluxes to 
their cata log values shows gene rally good agreem e nt. A s 



argued in VoUmer et ah] ( 2005 ) and Jacobs et al. (2011), 
there is considerable diifticulty in accurately comparmg 
catalogs, particularly at these frequencies. Differences in 
interferometer resolution and synthesized beam patterns 
can lead to non-trivial disagreement in measured source 
flux densities. Given this fact, and the complicated side- 
lobes associated with a 12-element array present in our 
method, we do not find the lack of better agreement with 
catalog fluxes troubling. 



To test the stability of our pipeline, we test a separate 
observation spanning the following day. The beam model 
produced by this data matches our first model to within 
2.5%, and source flux densities remain consistent to 5%. 
This test confirms our suspicion that our errors are not 
dominated by random noise. 

6. CONCLUSIONS 

We have presented a new technique for calibrating the 
primary beam of a wide-field, drift-scanning antenna ele- 
ment. The key inputs to the method are measurements of 
perceived flux density of individual sources versus time 
as they drift through t he beam. In this paper, we use 
delay/delay-rate fllters (Parsons & Backer 2009) as es- 
timators in a self-calibration loop to extract such mea- 
surements from raw interferometric data. However, the 
remainder of our method is agnostic to how these tracks 
are derived. The only assumption necessary to make this 
an over-constrained and solveable problem is 180° rota- 
tional symmetry in the beam. With this assumption, we 
create "crossing points" where there is enough informa- 
tion for a least-squares inversion to solve for both the 
inherent flux density of each source and the response of 
primary beam at the crossing points. 

We test this approach using simulated tracks of per- 
ceived source flux density across the sky and simulated 
visibilities. Using the source tracks, the least-squares 
inversion reliably recovers the primary beam values and 
source flux densities in the presence of noise 10 times that 
present in a 12-element PAPER array. The simulated 
visibilities demonstrate that sidelobes of other bright 
sources are the source of the dominant errors in source 
extraction when only 12 antennas are used. The presence 
of these features in the perceived source flux densities 
limits the accuracy of estimates of inherent source flux 
densities. However, in simulation, we are able to recover 
a primary beam model accurate to within 15% percent 
and source flux densities with an average error of 10%. 
Using prior information regarding beam smoothness, we 
improve our model to better than 10% accuracy. 

While these caveats about the effectiveness of the least- 
squares technique in the presence of sidelobes may seem 
worrisome, it bears repeating that these are issues with 
data quality and not with the technique itself. For exam- 
ple, data from a 32-element PAPER array has shown that 
DDR filters can extract the sources used in this analysis 
with little systematic biases. (We do not use this data 
here do the lack of a temperature record for gain sta- 
bilization and the presence of a particularly active Sun 
when the array was operating.) Therefore, it seems that 
this technique has significant potential for precise beam 
calibration on larger arrays. 

Another future goal for this technique is to calibrate 
the frequency-dependence of the beam. Here, we have 
used our entire bandwidth to improve signal-to-noise in 
source extractions. For a larger array with higher SNR 
and lower sidelobes, our perceived source flux density 
measurements can be cut into sub-bands to look at the 
beam as a function of frequency. 

Finally, it is possible to forgo the assumption of 180° 
rotational symmetry altogether and allow for possible 
north/south variation in the beam. An experiment in 
which the dipoles are ■physically rotated on a daily basis 
can be used to create the same kind of crossing points. 
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since one has changed the section of the beam each source 
crosses through. Work is progressing on such an experi- 
ment using the PAPER array. 
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Table 1 



Measured flux densities for all the sources used in this work. Catalog values come from the 3C catalog l |Edge et al.|1959| at 159 MHz 
unless otherwise noted. Errors arc estimated by measuring the change in flux required to increase A.x^ between the model and measured 

perceived flux density source tracks by 1.0. 



RA 


Dec 


Measured Flux Density (Jy) 


Catalog Flux Density (Jy) 


Name 


1:08:54.37 


+13:19:28.8 


78.1^24 4 


58, 49" 


3C 33 


1:36:19.69 


+20:58:54.8 


47.3;j?i 


27 


3C 47 


1:37:22.97 


+33:09:10.4 


6o.oti°;? 


50 


3C 48 


1:57:25.31 


+28:53:10.6 


oa 0+9.8 
■-^^ — 23 6 


7.5, 16\ 23= 


3C 55 


3:19:41.25 


+41:30:38.7 




50 


3C 84 


4:18:02.81 


+38:00:58.6 


75.9^2^io 


60 


3C 111 


4:37:01.87 


+29:44:30.8 


205.3j:^6-9 


204 


3C 123 


5:04:48.28 


+38:06:39.8 




85 


3C 134 


5:42:50.23 


+49:53:49.1 


56.5111, 


63 


3C 147 


8:13:17.32 


+48:14:20.5 




66 


3C 196 


9:21:18.65 


+45:41:07.2 


48.51^:^ 


42 


3C 219 


10:01:31.41 


+28:48:04.0 


40.3t»5^6 


30 


3C 234 


11:14:38.91 


+40:37:12.7 


M.5ll-9% 


21.5 


3C 254 


12:30:49.40 


+12:23:28.0 


1108.9t29-6 


1100 


Vir A 


14:11:21.08 


+52:07:34.8 


64.41^:^, 
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3C 295 


15:04:55.31 


+26:01:38.9 


73 5+^-^ 


72 


3C 310 


16:28:35.62 


+39:32:51.3 


60.5tli4 


49 


3C 338 


16:51:05.63 


+5:00:17.4 


461 7+^^-^ 


325^ 378\ 373^= 


Her A 


17:20:37.50 


-0:58:11.6 


228.8t?i7^5 


180, 236'^, 276^ 215^= 


3C 353 


18:56:36.10 


+1:20:34.8 


OKI O+60.0 


680 


3C 392^* 


19:59:28.30 


+40:44:02.0 


10622.9^}i■° 


10623 


Cyg A'^ 


20:19:55.31 


+29:44:30.8 


62 2+1^" 


36 


3C 410 


21:55:53.91 


+37:55:17.9 


^^■^-16. 2 


43 


3C 438 


22:45:49.22 


+39:38:39.8 


7-1 + 12.4 
' -'^■'-'-19.3 


50 


3C 452 


23:23:27.94 


+58:48:42.4 


9198 1+203.4 


6230^ 


Cas A 



" 3CRR dLaing ct al. U983[ l at 178 MHz 

Culgoora (Sice 199 ^^506 MHz. 
" PSA32 (Jacobs etg.|2011^ . Measurements made by PAPER 

Extended (35x27 arcminj supernova remn ant (W44) 



Flux density calibrator, using values from Baars et al. 
Ex trapolated using the fo rmula in [Baars et 
by |Helmboldt fc Kassim| l |2009^ . 



(19771. 

his IS probably a substantial underestimate of the flux density, as suggested 



